Preprint typeset in JHEP style - HYPER VERSION 



ITP-Budapest 568 
DESY 01-057 



Lattice determination of the critical point of QCD at 
finite T and fi 



Z. Fodor 

Deutsches Elektronen- Synchrotron DESY, Notkestr. 85, D-22607, Hamburg, Germany 
Institute for Theoretical Physics, Eotvos University, Pdzmdny 1, H- 1117 Budapest, 
Hungary. 

Email: f odor@poe . elte .hi^ 



S.D. Katz 

Institute for Theoretical Physics, Eotvos University, Pdzmdny 1, H- 11 17 Budapest, 
Hungary. 



Email:k.&tz©bodzi . elte .hu 



Abstract: Based on universal arguments it is believed that there is a critical point (E) 
in QCD on the temperature (T) versus chemical potential (/x) plane, which is of extreme 
importance for heavy-ion experiments. Using finite size scaling and a recently proposed 
lattice method to study QCD at finite [i we determine the location of E in QCD with 
nf =2+1 dynamical staggered quarks with semi-realistic masses on Lt = 4 lattices. Our 
result is Tg = 160 + 3.5 MeV and \xe = 725 + 35 MeV. For the critical temperature at 
H = we obtained T c = 172 ± 3 MeV. 
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1. Introduction. 

QCD at finite T and/or fj, is of fundamental importance, since it describes relevant features 
of particle physics in the early universe, in neutron stars and in heavy ion collisions (for a 
clear introduction see QCD is asymptotically free, thus its high T and high density 
phases are dominated by partons (quarks and gluons) as degrees of freedom rather than 
hadrons. In this quark-gluon plasma (QGP) phase the symmetries of QCD are restored. 
In addition, recently a particularly interesting, rich phase structure has been conjectured 
for QCD at finite T and /x % §. 

Extensive experimental work has been done with heavy ion collisions at CERN and 
Brookhaven to explore the [i-T phase diagram. Note, that past, present and future heavy 
ion experiments with always higher and higher energies produce states closer and closer to 
the T axis of the \i-T diagram. It is a long-standing open question, whether a critical point 
exists on the \i-T plane, and particularly how to predict theoretically its location H 

Let us discuss first the fi=0 case. Universal arguments f|| and lattice simulations 
§ indicate that in a hypothetical QCD with a strange (s) quark mass (m s ) as small as 
the up (u) and down (d) quark masses {m u ^) there would be a first order finite T phase 
transition. The other extreme case (nj=2) with light u/d quarks but with an infinitely large 
m s there would be no phase transition only an analytical crossover. Note, that observables 
change rapidly during a crossover, but no singularities appear (we will use the expression 
"transition" if we do not want to specify whether we deal with a phase transition or a 
crossover). This means that between the two extremes there is a critical strange mass {m c s ) 
at which one has a second order finite T phase transition. Staggered lattice results on L^=4 
lattices with two light quarks and m s around the transition T (nj=2+l) indicated |7| that 
m c s is about half of the physical m s . Thus, in the real world we probably have a crossover. 
(Clearly, more work is needed to approach the chiral and continuum limits. Note, that the 
puzzle due to an unexpected strengthening observed § for the nj=2 Wilson action with 
intermediate m U ( i was resolved by using an improved action ||].) 
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Returning to a non- vanishing /i, one realizes that arguments based on a variety of 
models (see e.g. [[R], [2|, [|) predict a first order finite T phase transition at large \x. Com- 
bining the /i = and large /i informations an interesting picture emerges on the \i-T plane. 
For the physical m s the first order phase transitions at large [i should be connected with 
the crossover on the [i = axis. This suggests that the phase diagram features a critical 
endpoint E (with chemical potential [ie and temperature Te), at which the line of first 
order phase transitions (fi > fiE and T < Te) ends |B|. At this point the phase transition 
is of second order and long wavelength fluctuations appear, which results in character- 
istic experimental consequences, similar to critical opalescence. Passing close enough to 
(/J:E,Te) one expects simultaneous appearance of signatures (e.g. freeze-out type behavior 
of observables constructed from the multiplicity and transverse momenta of charged pions), 
which exhibit nonmonotonic dependence on the control parameters [11], since one can miss 
the critical point on either of two sides. 

The location of this critical point is an unambiguous, non-perturbative prediction of 
the QCD Lagrangian. Unfortunately, no ab initio, lattice analysis based on QCD was done 
to locate the endpoint. Crude models with m s = oo were used (e.g. ||) suggesting that 
[iE ~ 700 MeV, which should be smaller for finite m s . The result is sensitive to m s , thus 
for realistic cases previous works could not predict the value of \xe even to within a factor 
of 2-3. The goal of this exploratory work is to show how to locate the endpoint by a lattice 
QCD calculation. We use full QCD with dynamical nf =2+1 staggered quarks. 

QCD at finite \x can be formulated on the lattice [12]; however, standard Monte-Carlo 
techniques can not be used at fi ^ 0. The reason is that for non-vanishing real [i the func- 
tional measure -thus, the determinant of the Euclidean Dirac operator- is complex. This 
fact spoils any Monte-Carlo technique based on importance sampling. Several suggestions 
were studied to solve the problem. Unfortunately, none of them was able to locate (/ie,Te)- 

In a recent paper we proposed a new method ]13] to study lattice QCD at finite T and 
[i. The idea was to produce an ensemble of QCD configurations at /x=0 and at T c . Then we 
determined the Boltzmann weights ( |14fl of these configurations at fi ^ and at T lowered 
to the transition temperatures at this non-vanishing /i. Since transition configurations 
were reweighted to transition configurations a much better overlap was observed than 
by reweighting pure hadronic configurations to transition ones |0|. We illustrated the 
applicability of the method in nj=4 dynamical QCD. Using only O(10 3 -10 4 ) configurations 
quite large fi could be reached and the transition line separating the hadronic phase and 
the QGP was given on the \x-T plane. 

In this letter we generalize the above method to arbitrary number of staggered quarks. 
We apply it to the nf =2+1 case. We determine the volume (V) dependence of the Lee- Yang 
zeros of the partition function on the complex gauge coupling {(3) plane. Based on this V 
dependence we determine the type of the transition as a function of fi. The endpoint {j,e 
is given by the value at which the crossover disappears and finite-V scaling predicts a first 
order phase transition. These finite T calculations are done on Lt = 4 lattices. In order to 
set the physical scale we determine the pion and rho masses (m^, m,,), the string-tension 
(a) and the Sommer [[Hj] scale (Rq) at T = 0. Our quark masses are "semi-realistic": m s is 
set about to its physical value, whereas m u a are approximately four times as heavy as they 
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are in the real world. Having determined the lattice spacing we transform our result to 
physical units and give T c the location of (/j,e,Te) and show the phase diagram separating 
the hadronic phase and the QGP. 

Though this study performes a V— > oo extrapolation, larger volumes, larger Lj-s and 
smaller masses are also needed to give the final answer to (he^Te)- 



2. Staggered quarks at jj, ^ 0. 



The partition function of QCD with nj degenerate staggered quarks (for an introduction 
see eg. [p] ]) is given by the functional integral of the bosonic action S b at gauge coupling 
f3 over the link variables U, weighted by the determinant of the quark matrix M, which 
can be rewritten pT 



as 



Z{p,m,n) = J VUexp[-Sb(l3,U)}{detM{m,n,U)} nf/i 
VU exp[-S b ((3 w , U)\ [det M(m, » w , U)] n f' A 

[det M(m, fi, U)] n f/ A \ 



(2.1) 



exp[-S b (0,U) + S b (p w ,U)] 



[det M(m, Hw, U)] n f/ 4 



where m is the quark mass, \x is the chemical potential of the quark. For non-degenerate 
masses one uses simply the product of several quark matrix determinants on the 1/4-th 
power. Standard importance sampling works and can be used to collect an ensemble of 
configurations at (3 W and \x w with Re(/i to )=0. It means we treat the terms in the curly 
bracket as an observable -which is measured on each independent configuration- and the 
rest as the measure. By simultaneously changing j3 and /j, one can ensure that even the 
mismatched measure at (3 W and [i w samples the regions where the original integrand with 
(3 and n is large. In practice the determinant is evaluated at some \x and a Ferrenberg- 
Swendsen reweighting ]l4j is performed for the gauge coupling (3. 

Due to the complex nature of det M(m, /i, U) an additional problem arises, one should 



choose among the possible Riemann-leaves of the fractional power in eq. (2.1). This can be 
done by using the fact that at fj, = fj, w the ratio of the determinants is 1 and the ratio should 
be a continuous function of fi. However, the continuity can only be ensured if the analytical 
dependence of the determinant on fi is known (the determinant oscillates strongly with fj,, 
so measuring it for several \x values is not satisfactory). This dependence can be given by 
the following way (the idea goes back to a method of Ji~8|]). 
First gauge fix to Aq = on all but the last timeslice 



/ #0 







M 



-e-» Bo 



(2.2) 



\-ePT 
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The Bi are 3-Ljj x 3Lf matrices containing the mass and spatial hopping terms (3 is the 
number of colors) and T contains the only remaining temporal links on the last timeslice. 
By multiplying the j-th column of M by e~ m and the j-th row by e- 7 ^ and moving the 
leftmost column to the right one gets a matrix with the same determinant 



/ 1 o 

Bi 1 



V 



1 

B L t -l 



Bo 
-1 



-e Ltfi T 



\ 



(2.3) 



We evaluate the determinant by Gauss-elimination. After Lt — 2 steps we get a 6Ljj x 6Ljj 
matrix 

/ 1 + ci • e~ Ltfl c 2 \ 

1 c 3 + c 4 • e~ Lt ^ c 5 - e Lt ^T I ' 



(2.4) 



where q are //-independent matrices. It is straightforward to show that 

det M = e - 3i ' Lt ^ det(P - e Ltl "), 



where 



C3C1 - c 4 (c 5 - c 3 c 2 )Tt 



(2.5) 



(2.6) 



This is just the characteristic equation of P. To find the Aj eigenvalues one needs 0{L S ) 
operations, which gives the determinant as a function of fi and determines the ratio of the 
fractional powers in eq. ( |2.1[) unambiguously 



det M(n) = e 



-3LU 



6L| 

IT 



(2.7) 



Using eq. ( p.l| ) with the determinant given by eq. (2/7) we can give the partition 
function for complex /i and (3 values. In the following we keep [i real and look for the zeros 



of the partition function on the complex (3 plane. These are the Lee- Yang zeros [19], whose 
V ^ 00 behavior tells the difference between a crossover and a first order phase transition. 
At a first order phase transition the free energy oc log Z((3) is non-analytic. Clearly, a 
phase transition can appear only in the 00 limit, but not in a finite V. Nevertheless, 
the partition function has zeros at finite V, the Lee- Yang zeros, at "unphysical" complex 
values of the parameters, in our case at complex (3. For a system with a first order phase 
transition these zeros approach the real axis in the V— > 00 limit (detailed analyzes suggest 
1/V scaling). This 00 limit generates the non-analyticity of the free energy. For a 
system with crossover the free energy is analytic, thus the zeros do not approach the real 
axis in the V— * 00 limit. The Lee- Yang technique was successfully applied to determine 
the endpoint of the electroweak phase transition (2(], 21]. 
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3. T ^ and T = results for n f =2+l. 



Using the formulation described above we study nf =2+1 QCD at T ^ on Lj = 4, L s = 
4, 6, 8 lattices and at T = on V = 10 3 • 16 lattices. m u ^ = 0.025 and m s = 0.2 were 
chosen for the bare quark masses. We used the R algorithm of the MILC collaboration's 



code |22]. 





0.1 


0.2 


0.3 


0.4 


Re(Po); L s = 4 
10 2 Im(/? ) 


5.151(1) 
5.56(8) 


5.141(1) 
5.50(9) 


5.127(2) 
5.42(15) 


5.121(5) 
5.56(38) 


Re(/3o); L s = 6 
10 2 Im(/? ) 


5.193(1) 
2.66(6) 


5.174(1) 
2.54(9) 


5.152(3) 
2.19(31) 


5.143(7) 
1.82(39) 


Re(/? ); L s = 8 
10 2 Im(/? ) 


5.193(1) 
1.38(6) 


5.172(1) 
1.32(17) 


5.159(1) 
1.31(18) 


5.140(1) 
0.48(7) 


Re(/? );^s -> oo 
10 2 Im(/? ) 


5.201(1) 
1.02(6) 


5.178(1) 
1.12(11) 


5.162(2) 
0.74(19) 


5.143(2) 
-0.25(10) 


/3 


m n 


m p 


Rq 




5.208 
5.164 


0.393(2) 
0.393(2) 


1.22(2) 
1.28(3) 


1.87(3) 
1.76(5) 


0.58(7) 
0.75(5) 



Table 1: T ^ and T = results. The upper part is a Summary of the Lee- Yang zeros obtained 
at different chemical potentials, while the lower part shows the measured T = observables for two 
(3 values. 
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At T 7^ we determined the com- 
plex valued Lee- Yang zeros, (3q, for dif- 
ferent V-s as a function of fi. Their 
oo limit was given by a Po(V) = + 
C/V extrapolation. The results (listed in 
Table 0) are from 14000, 3600 and 840 
configurations on L s =4,6 and 8 lattices, 
respectively. Figure shows Im(/?g°) 
a function of fi enlarged around the end- 
point Hend- The picture is simple and 
reflects the physical expectations. For 
small fi-s the extrapolated Im(/3^°) is in- 
consistent with a vanishing value, and 
the prediction is a crossover. Increasing fi the value of Im(/?Q°) decreases, thus the tran- 
sition becomes consistent with a first order phase transition. (Note, that errors decrease 
close to the endpoint, and the Im(/?Q°) extrapolation, due to the relatively small volumes, 
slightly overshoots. Both phenomena were observed already in the electroweak case e.g. 
1 21 ] ) . The statistical error was determined by a jackknife analysis using subsamples of the 
total L s = 4, 6 and 8 partition functions. The systematic uncertainty, estimated from the 
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Figure 1: 

potential. 



Im(/3o°) as a function of the chemical 
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overshooting, was added linearly to the statistical error. The dashed line of Figure |l| shows 
the fit and leads to our primary result; l^end — 

0.375(20). 

Table |] contains also the T = results. To set the physical scale we used a weighted 
average of R (1/403 MeV), m p (770 MeV) and yfa (440 MeV), obtained in lattice units 
by different fitting procedures. It is important to note, that (including systematic errors 
due to finite V) we have (Ro ■ m n ) = 0.73(6), which is at least twice as large as its physical 
value. Thus our m u ^ is at least four times larger than it should be. 

Let us estimate the applicability of 
the method approaching the chiral and 
continuum limits. In the present analy- 
sis the evaluation of the eigenvalues was 
somewhat less costly than the produc- 
tion of the configurations. For physical 
m u 4 the latter would need an additional 
factor of O(50) (the former remains the 
same). Thus, for physical masses at least 
upto V=4-12 3 the cost of the eigenvalue 
determination is subdominant. Extend- 
ing the analysis to this volume reduces 
the error on fi en d to a level, which is 
not even needed (uncertainties due to fi- 
nite lattice spacing could be more impor- 
tant). Since for finer lattices the eigen- 
value evaluation goes with and the configuration production at least with Lj? the eigen- 
value evaluation remains subdominant. At physical masses fi e nd is probably closer to the 
fi=0 axis (for recent lattice works see M). Thus, the overlap between /i=0 and /i / 
configurations is even better. It means less statistic might be enough to apply eq. (|2.1| ) 
than it was used in this work. Note, that the quark masses of the present work are half of 
those used in Ref. [jnj; however, in both cases it was possible to reweight in /i far beyond 
m n /2 (the typical premature onset /i value of the Glasgow method |T5[] ). Thus, we expect 
that our method does not suffer from this type of onset problem when approaching the 
chiral limit. 

Figure || shows the phase diagram in physical units, thus T as a function of fj,B, the 
baryonic chemical potential (which is three times larger then the quark chemical potential). 
The endpoint is at T E = 160 ± 3.5 MeV, \k E = 725 ± 35 MeV. At fi B =0 we got T c = 
172 ± 3 MeV. 

4. Conclusions, outlook. 

We used a recently proposed method |l3|] and studied the \x-T phase diagram of QCD 
with dynamical nj =2+1 quarks. We presented an ah initio technique to determine the 
location of the endpoint. Using the above method we obtained Tg=160±3.5 MeV for the 
temperature and /i£=725±35 MeV for the baryonic chemical potential of the endpoint. 
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Figure 2: The phase diagram in physical units. 
Direct results are with errorbars. Dotted line illus- 
trates the crossover, solid line the first order phase 
transition. The small box shows the uncertainties of 
the endpoint. 
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This result was based on the oo behavior of the Lee- Yang zeros of the partition 
function. We used L t =A and our quark masses were "semi-realistic" (m s was set to about 
its physical value, whereas m U) d were four times heavier than in the real world). Though 
He is too large to be studied at RHIC/LHC, the endpoint would probably move closer to 
the //=0 axis when the quark masses get reduced. At fi=0 we obtained T c =172±3 MeV. 
Clearly, more work is needed to get the final values. One has to extrapolate to zero step-size 
in the R-algorithm and to the thermodynamic, chiral and continuum limits. 

We thank F. Csikor, F. Karsch, I. Montvay and A. Ukawa for useful comments on the 
manuscript. This work was partially supported by Hungarian Science Foundation grants 
No. OTKA-34980/29803/22929/28413/OMMU-708/IKTA. This work was in part based 



on the MILC collaboration's public lattice gauge theory code [22]. 
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